
function [ES,S_2S,DES] = surplus_2S(gamma,B,B2S,D,P,xi,n,M)

N0 = n(1,1); N1 = n(1,2); N2 = n(1,3); N = N0 + N1 + N2;
NS = size(xi,1);
M = M - 1; % party affiliation of coalition candidate (for menu fixed effects)
need_grad = nargout > 2;

% PRI-PVEM joint surplus from voting second stage (and gradient)
ES = zeros(N,NS);
S_2S = zeros(N,3,NS);
if need_grad
    DES = zeros(N,length(gamma)+length(B)+length(B2S),NS);
end

for t = 1:N
    Xt = [kron(eye(2),[1-M,M,D(t,:)]),log(P(t,3:4)')];
    parfor ns = 1:NS
        es = Xt * B2S + xi(ns,:)';
        es = exp(es);
        es = es ./ (1 + sum(es,1));
        ES(t,ns) = gamma(3) * log(es(1,:) + 0.5 * (1 - sum(es,1))) + gamma(4) * log(es(2,:) + 0.5 * (1 - sum(es,1))); %#ok<PFBNS> 
        St(:,ns) = [es(1,:);es(2,:);1-sum(es,1)];
        if need_grad
            DEStns = zeros(length(gamma)+length(B)+length(B2S),1);
            DEStns(3) = log(es(1,:) + 0.5 * (1 - sum(es,1)));
            DEStns(4) = log(es(2,:) + 0.5 * (1 - sum(es,1)));
            for k = 1:size(Xt,2)
                xx = Xt(:,k);
                sx = sum(es .* xx,1);
                sx = [es(1,:) .* (xx(1,:) - sx);es(2,:) .* (xx(2,:) - sx)];
                DEStns(length(gamma)+length(B)+k) = ...
                    gamma(3) * ((sx(1,:) - 0.5 * sum(sx,1)) ./ (es(1,:) + 0.5 * (1 - sum(es,1)))) +...
                    gamma(4) * ((sx(2,:) - 0.5 * sum(sx,1)) ./ (es(2,:) + 0.5 * (1 - sum(es,1))));
            end
            DESt(:,ns) = DEStns;
        end
    end
    S_2S(t,:,:) = St;
    if need_grad
        DES(t,:,:) = DESt;
    end
end



